!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine model_variables (ro, mx, my, mz, bx, by, bz, en, gx, gy, po, x, y, nx, ny)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny

   real(8), dimension(nx), intent(in) :: x
   real(8), dimension(ny), intent(in) :: y
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy, po

   integer(4) :: i, j

   real(8) :: current_sheet_center_x
   real(8) :: current_sheet_center_y
   real(8) :: current_sheet_len_x
   real(8) :: current_sheet_len_y
   real(8) :: pr

!----------------------------------------------------------------------------------------------------------------------------------|
!  Current sheet
!----------------------------------------------------------------------------------------------------------------------------------|
   current_sheet_center_x = 0.0d0
   current_sheet_center_y = 0.0d0
   current_sheet_len_x = 0.005 * rangex
   current_sheet_len_y = 0.005 * rangey

   do j = 1, ny
   do i = 1, nx

      pr = 1.0d0
      ro(i, j) = gm
      mx(i, j) = 0.0d0
      my(i, j) = 0.0d0
      mz(i, j) = 0.0d0
      by(i, j) = 0.0d0
      if (y(j) - current_sheet_center_y .lt. -current_sheet_len_y) then
         bx(i, j) = -1.0d0
         bz(i, j) = 0.0d0
      elseif (y(j) - current_sheet_center_y .gt. current_sheet_len_y) then
         bx(i, j) = 1.0d0
         bz(i, j) = 0.0d0
      else
         bx(i, j) = 1.0d0 * sin ((y(j) - current_sheet_center_y) / current_sheet_len_y * pi / 2.0d0)
         bz(i, j) = 1.0d0 * cos ((y(j) - current_sheet_center_y) / current_sheet_len_y * pi / 2.0d0)
         mx(i, j) = 0.05d0 * ro(i, j) * sin (x(i) / rangey * 16.0d0 * pi) * sin (y(j) / current_sheet_len_y * pi)
      endif
      en(i, j) = (bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)) / 2.0d0 + pr / gmm1
      gx(i, j) = 0.0d0
      gy(i, j) = 0.0d0

      po(i, j) = 0.0d0

   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
end subroutine model_variables
